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SUMMARY 


Problems requiring the determination of eigenvalues occur frequently in the study 
of dynamical systems. The formulation and processing of the equations describing 
these systems involve time-consuming matrix operations. Although numerical methods 
are available, it is more convenient and informative to formulate the relevant 
matrices and transfer functions symbolically. This can be achieved by the use of 
symbolic computation, which is a technique whereby algebraic operations, symbolic dif- 
ferentiation, matrix formulation and inversion, etc., can be performed on a digital 
computer equipped with a formula— manipulation compiler. An example is included that 
demonstrates the facility with which the system dynamics matrix and the control dis- 
tribution matrix from the state space formulation of the equations of motion can be 
processed to obtain eigenvalue loci as a function of a system parameter. 


INTRODUCTION 


The response of a dynamical system to control inputs may be stable or unstable, 
depending on the characteristics of the system parameters. An Important aspect of 
stability studies is the determination of the eigenvalues of the system dynamics 
matrix in the case of a continuous system, and the eigenvalues of the corresponding 
transition matrix in the case of discrete systems. It is shown that symbolic compu- 
tation facilitates the determination of the eigenvalues of open-loop systems and the 
data required to plot stability loci for closed— loop systems. 

The locations of the eigenvalues in the complex plane provide the researcher with 
a clear indication of the stability of the system being considered and the extent to 
which it is damped. Moreover, an examination of the eigenvalue loci suggests what 
steps must be taken to ensure a satisfactory response. 

The example chosen to demonstrate the technique is a fourth-order system repre- 
senting the longitudinal response of -a DC 8 aircraft to elevator inputs. This simpli- 
fied system has two dominant modes, one of which is lightly damped and the other well 
damped. The loci may be used to determine the value of the controlling parameter 
that satisfies design requirements. 

The results were obtained using the MACSYMA Symbolic Manipulation System. This 
is a large computer-programming system used for performing symbolic as well as numer- 
ical mathematical manipulations. It was developed by the Mathlab Group of the MIT 
Laboratory for Computer Science. 



ANALYSIS 


Laplacfi Transforms and Transfer Functions 

State space equations- For a linear, time -invariant dynamical system, which is 
free from disturbance inputs, the vector-matrix form of the state space equations is 
(ref. 1) 

X = AX + BU ' (1) 

where X is the state vector and U is the input vector. In general, A and B are 
nxn and nxm matrices, respectively. If the output vector Y is a linear combina- 
tion of the state vector and the input vector, it may be written as follows; 

Y = CX + DU (2) 

Here, C- and D are pxn and pxm constant matrices, respectively. The general system 
described by equations (1) and (2) has n state variables, m inputs, and p outputs. 
See figure 1. 

STATE 

input variables 



Figure 1.- The state variable formulation for a multi-input, multi-output system. 

From equation (1), with zero initial conditions 

. / (si - A)X(s) = BU(s) (3) 

where X(s) and U(s) are the Laplace transforms of the state vector and the input 
vector, respectively. The Laplace transform of the state vector is obtained by solv- 
ing equation (3) 

X(s) = (si - A)-^BU(s) (A) 


From equation (2) 


Y(s), = CX(s) + DU(s) 


By assuming that D is a null matrix the output simplifies accordingly, and 
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Y(s) = CX(s) 


(5) 


Determination of Eigenvalues in a Practical Case 

Onen-loop transfer functions- A state space formulation of the equation of motion 
describing the longitudinal response of a DC 8 aircraft to an elevator input, has the 
following form (ref. 2). 
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where U is the elevator displacement. In this case, the output Y is 


Y = c"^X 


= (Bo Bi 1 0) 


Xi 

X2 

X3 

LX4, 


(7) 


Before transfer functions for the individual states can be obtained, equation (4) 
requires that the identity matrix, the Laplace operator s, and the system dynamics 
matrix A, be combined as shown. This operation is described in terms of reproduced 
computer output as follows; 

First, the system dynamics matrix A is entered 

r 0 1 0 0 1 


I 0 0 0 1 

L-Aq -Ai -A2 -A3J 

Next, the computer provides the fourth-order identity matrix I and multiplies 
it by s 
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When equation (8) is subtracted from equation (9) the matrix (si - A) is 
obtained 



( 10 ) 


Laq Aj A2 s + A3 


Inversion of this matrix yields the coefficient of BU(s) in equation (4). 
Because of space limitation, the result is printed out in four columns (ref. 3). 


Col 1 = 


Col 2 = 


$3 + A3S2 + A2S + Ai 
s'* + A3S3 + A2S2 + AiS + Aq 

^0 

s'* + AgS^ + A2S2 + AjS + Ag 

AqS 

s'* + AgS 3 + A2S2 + AjS + Aq 

AqS2 

s'* + + A2S2 + AjS + Aq 

S2+A3S + A2 

s'* + AjS^ + A2S^ + AjS + Ag 
■ S^ + A3S^ + A 2 S 
s'* + A3S3 + A2S2 + AiS + Aq 
A jS + Aq 

s'* + AsS^ + A2S2 + AiS + Aq 



( 11 ) 


A^S^ + AgS 

s'* + AjS^ + A2S2 + AjS + Aq 

i 
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Col 3 


S + A3 


s'* + AaS^ + A2S2 + AjS + Ao 


$2 + A3S 


s'* + A3S3 + A2S2 + AjS + Aq 


+ A3S2 


s'* + AgS^ + A2S2 + AjS + Aq 


A2S2 + AjS + Aq 


s'* + AsS^ + A2S2 + AiS + Ao 



Col 4 = 


1 


s'* + AjS^ + A2S2 + AjS + Aq 


S 


s'* + A3S^ + A2S2 + AjS + Aq 


S 2 


s'* + AjS^ + AgS^ + AjS + Aq 


S^ 


s'* + AgS^ + A2S2 + AjS + Aq 


- 



( 11 ) 

Concluded 


Before equation (4) can be expanded, the control distribution matrix must be 
entered. 


0 

0 

0 

K 


= B 


( 12 ) 


Premultiplication of this vector by (si - A)“^ yields the following column vector; 
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K 


+ AgS^ + A2S2 + A^S + Aq 

, KS 

s'* + A'aS^ + A2S2 + AiS + Aq 

■ KS^ 

s'* + A3S3 + A2S2 + AiS + Aq 

KS3 

s'* + AgS^ + A2S2 + A^S + Aq 


( 13 ) 


If zero initial conditions are assumed, the Laplace transform of the input is U(s) . 


Multiplication of equation (13) by U(s) gives the Laplace transforms of the 
states. That is: 


X(s) 


KU(S) 

s'* + A3S3 + A2S2 + AjS + Aq 

KSU(S) ^ 

s'* + AjS^ + A2S2 + AjS + Aq 

KS^UCS) 

s'* + AgS 3 + A2S2 + AjS + Aq 

KS^U(S) 

s'* + AgS^ + A2S2 + AjS + Aq 


(14) 


Up to this point the only operations undertaken were addition, subtraction, mul- 
tiplication, and inversion of matrices. These operations require the following sim- 
ple program statements. 

Addition and subtraction of the matrices A and B: A ± B; 

Multiplication: A.B; 

Inversion of the matrix A: AAA-1; 

Two programming steps are required to obtain the Laplace transforms of the indi- 
vidual states. These are: 

FOR 1:1 THRU 4 DO ROW[I] :FIRST (ROW( (D7) ,1) )$ 

FOR 1:1 THRU 4 DO (X[I] :R0W[1][1], DISPLAY (X[l])); 

The resulting Laplace transforms Xj^ are 
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^ KU(S) 

^ + A2S^ + AiS + Ao 

^ ^ KSU(S) 

^ s'* + A3S3 + A2S2 + AiS + Aq 

^ KS^U(S) 

^ s*t + AgS^ + A2S2 + A^s + Aq 

^ KS^UCS) 

~ s'* + AgS^ + A 2 S 2 + AjS + Aq 


\ 




/ 


( 15 ) 


The corresponding transfer functions for the single input case are obtained 

by using the following program statements; 

FOR I;1 THRU 4 DO T[I]: (X[I]/U(S))$ 

FOR 1:1 THRU 4 DO DISPLAY(T[I]) ; 


Ti = 


K 


^ s'* + AoS^ + A 2 S 2 + AjS + A^ 


To = 


KS 


s'* + AoS^ + A,S^ + A,S + Af 


T3 = 


Tu = 


KS" 


s'* + A3S3 + A2S2 + AjS + Aq 
KS^ 

s'* + AjS^ + + AjS + A(. 


(16) 


The Laplace transform of the output Y(s) is obtained by expanding equation (7). 


where 


Y(s) = C^X(s) 

= (Bo Bi 1 0) 


That is 

Y(s) = c'*^(sl - A)“^BU(s) 

Substitution of the vector (14) in this equation gives the Laplace transform of the 
output 


Y(s) 


(KS^ + BiKS + BqK)U(S) 

s'* + AgS^ + A2S2 + A^S + Aq 
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The corresponding transfer function T(s) is 

' KS 2 + B,KS + Bf,K 

X(s) = ^ ^ =— ^ 

s'* + A3S3 + A2S2 + A^S + Aq 

Open-loop eigenvalues- The eigenvalues of the open-loop system are the roots of 
the determinant of the matrix (si - A), or alternatively, the poles of the open-loop 
transfer function. Subsequent to the formulation of (si - A), the determinant is 
obtained as follows; 

DETERMINAlirr(%) ; 


S(S(S(S + A3) + A2) + Aj) + Aq 
The expanded form of this expression is obtained by using 
EXPAND(%); 

s'* + A 3 S 3 + A 2 S 2 + AjS + Aq 

The eigenvalues of the open-loop system are the roots of this equation. These 
can be evaluated by using the root finding program. However, before the root finding 
program can be implemented, all the coefficients must be replaced by their numerical 
values. This is easily accomplished by using a substitution routine. The following 
are values for the landing approach phase of the motion (also discussed in ref. 2). 


Bo 

;= 0.0323675 

Ao = 0.07006953 

Bl 

= 0.5955 

Ai = 0.09712526 



A 2 = 2.6813872 



A 3 = 1.700522 


Substitution of these values in the determinants! equation yields: 

s'*/*- 1.700522 s3 + 2.6813872s 2 + 0.09712526S + 0.07006953 
Using the notation %I to denote the complex number /^, the roots are 

S = 0.163195707 %I - 9.9571822E-3 
S = -0.163195707 %I - 9.9571822E-3 
S = 1.38386287 %I - 0.84030382 
S = -1.38386287 %I - 0.84030382 

These are the eigenvalues of the open-loop system. Since the ranges of the closed- 
loop eigenvalue loci are determined by the zeros of the transfer functions, these must 
be determined. The only transfer function of interest in the present application is 
T(s), the zeros of which are the roots of the equation 
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(s^ + 0.5955s + 0.0323675) = 0 


which are 


s = -0.0604999, s = -0.535000 

It should be noted that In addition to the two finite zeros, there are two zeros at 
infinity. 

Closed-toop transfer functions- In addition to the plant input, closed-loop 
systems have a reference or system input R(t). See figure 2. 


SYSTEM PLANT r 


PLANT 


PLANT 

AND 

SYSTEM 



Figure 2.- Single-input, single-output state 
variable feedback system. 


For this reason, equation (1) will be rewritten as follows: 

X = AX + B(tj + BiU + BqU) 


This formulation of the state equations requires that equation (7) be modified accord- 
ingly. That is 




Y = C% = (1 0 0 0) 


X2 

X3 



If initial conditions are assumed to be zero, the corresponding Laplace transform 
equations are 

sX(s) = AX(s) + B(s 2 + Bis + Bq)U(s) 
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therefore 


X(s) = (si - A)-1 b(s 2 + Bjs + Bq)U(s) 

The Laplace transform of the system input R(s) is related to the Laplace transforms 
of the plant input U(s) and the state X(s) by the equation 

U(s) = R(s) - kX(s) 

However, since the intention here is not to emphasize the control system aspect of the 
problem, but rather to demonstrate the feasibility of using symbolic computation to 
find the eigenvalues of dynamical systems, this equation will be simplified as 
follows: 

U(s) = R(s) - X(s) 

Substitution in the preceding equation yields the Laplace transform of the closed-loop 
state 

X(s) =1^1+ (s^ + BjS + BqHsI - A)~^bJ \s2 + BjS + Bq)(sI - A)-^BR(s) (17) 

where I is a (4 x 4) identity matrix. 

To permit inversion of the matrix 



+ Bj s + Bq) (sI 



(18) 


the (4 X 1) column vector B must be converted to a (4x4) matrix as follows: 


0 0 0 0 

0 0 0 0 

0 0 0 0 

K 0 0 0 


(19) 


With this modification, the required Laplace transforms are obtained. First, equa- 
tion (19) is multiplied by (s^ + Bjs + Bq) 


0 

0 

0 


K(s 2 + BiS + Bo) 


0 0 0 
0 0 0 
0 0 0 
0 0 0 


( 20 ) 


Premultiplication of matrix (20) by matrix (11) yields the following matrix: 
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K(s 2 + BiS + Bo) 

0 0 0 

+ A3S^ + A 2 S 2 + AjS + Aq 

KS(s 2 + BiS + Bo) 


s'* + A3S3 + A^S^ + AjS + Ag 

KS^(S^ + BiS + Bp) 

s'* + AgS^ + + AjS + Aq 


KS^(S^ + BiS + Bp) 
s'* + A 3 S 3 + A 2 S 2 + AjS + Aq 


( 21 ) 


When the identity matrix is added to matrix (21), the expanded form of matrix (18) is 
obtained. 


K(s2 + B^S + Bq) 
s'* + A3S3 + A2S2 + A^S + Aq 
KS(s 2 4. B^s + Bq) 
s'* + AgS^ + A2S2 + AjS + Aq 
KS 2 (S^ + BiS + Bo) 
s'* + A3S3 + A2S2 + AjS + Aq 
KS3(S^ + BiS + Bq) 
s'* + AgS^ + A2S2 + AjS + Aq 


0 0 0 
10 0 
0 10 
0 0 1 


( 22 ) 


Before equation (17) can be expanded, the matrix (22) must be inverted. The inverted 
form of this matrix is 


s'* + AaS^ + A2S^ + AiS + Aq 

s'* + AgS^ + (K + A2 )s 2 + (BjK + Aj)S + BqK + Aq 

KS^ + BiKS^ + BqKS 

s'* + A3S3 + (K + A2 )s 2 + (BiK + Ai)S + BpK + Ap 

KS ** + BjKS^ + BqKS 2 

s'* + AsS 3 + (K + A2)s 2 + (BiK + Ai)S + BpK + Ap 

KS^ + BiKS'* + BpKS^ 

s'* + AaS^ + (K + A2)S^ + (BiK + Ai)S + BqK + Ap 


0 0 
0 0 
1 0 
0 1 


(23) 


Premultiplication of matrix (21) by matrix (23) yields the matrix 
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KS^ + B^KS + BqK 

s'* + A3S3 + (K + A2 )s 2 + (BiK + Ai)S + BqK + Aq ° ° ° 

KS^ + BiKS^ + BqKS' 

S'* + AgSS + (K + A2 >s 2 + (BjK + Ai)S +, BqK + Aq ° ° ° 

KS** + B,KS^ + BqKS^ - 

s'* + AgS^ + (K + A2 >s 2 + (BjK + Aj)S + BqK + Aq ^ ^ ^ 

KS^ + BjKS'* + BqKS^ 

s'* + A3S3 + (K + A2)s2 + (BiK + Ai)S + BqK + Aq ° ° ° 


( 24 ) 


Formulation of the reference input vector which has only one component Ri (s) and pre- 
multiplication by matrix (24) completes the expansion of equation (17), and gives the 
column vector of closed-loop Laplace transforms for the individual states. 


Rl(S) 

0 

0 

0 

_J (KS^ + BjKS + BoK)Ri(S) 

s'* + A 3 S 3 + (K + A2)s 2 + (B]K + Ai)S + BqK + Aq 

(KS3 + BjKS^ + BqKS)Ri(S) 

s'* + AgS^ + (K + A2)s2 + (BjK + Aj)S + BqK + Ag 

(KS** + BiKS3 + BqKS^)R,(S) 

s'* + AgS^ + (K + A2)s2 + (BjK + Aj)S + BqK + Aq 

(KS^ + BiKS** + BqKs3)Ri(S) 

_S'* + A 3 S 3 + (K + A2)s2 + (BiK + Ai)S + BqK + Aq 


(25) 


(26) 


The program statements used to obtain the Laplace transforms of the individual open- 
loop states can be used in this case also. The components of the closed-loop vector 
are 
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( 27 ) 


(KS^ + B^KS + BqK)Ri (S) 1 

“ s'* + A 3 S 3 + (K + A2>s 2 + (BiK + Ai)S + BqK + Aq 

(KS3 + BiKS 2 + DgKS)JliCS> 

~ s'* + A 3 S 3 + (K + A2>s 2 + (BjK + Ai)S + BqK + Aq 
(KS** + BjKS^ + BqKS2)Rj^(S) 

" s'* + A 3 S 3 + (K + A2>s 2 + (BjK + Ai)S + BqK + Aq 

(KS5 + B^KS** -I- BqKS3)Ri(S) 

" s'* + A 3 S 3 + (K + A2>s 2 + (B^K + Ai)S + BqK + Aq 

closed-loop transfer functions, which are defined by (C39) , and the 
activating program statements are 

(C39) FOR 1:1 THRU 4 DO T[I]: (X[I]/R(S))$ 

(C40) FOR 1:1 THRU 4 DO DISPLAY(T[ I] ) ; 

KS2 + BjKS + BqK 

" s’** + A 3 S 3 + (K + A2>s2 + (BjK + Aj)S + .BqK + Aq 
KS3 + BjKS2 + BqKS 

^ — - ■ ^ ■ 

^ s'* + A 3 S 3 + (K + A 2 >s 2 + (BjK + Ai)S + BqK + Aq 

KS** + BiKS3 + BqKS2 

~ s'* + A 3 S 3 + (K + A2>s 2 + (BiK + Ai)S + BqK + Aq 

KS3 + BiKS** + BqKS3 

" s'* + A 3 S 3 + (K + A 2 >s 2 + (BjK + Aj)S + BqK + Aq 


Xl 

X2 

X3 

X4 

The corresponding 


The determinant of the matrix (22) is obtained by using the program statement 


DETERMINANT (22); 


K(s 2 + B^S + Bq) 


s'* + A3S3 + A2S2 + AjS + Aq 


+ 1 


The closed-loop eigenvalues are the roots of the determinant | (22) | - 0 


K(s 2 + BiS -t- Bp) 
s'* + A3S3 + A2S2 + AjS + Aq 


+ 1 = 0 
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that is, the roots of . 

I' ■ 

s'* + A 3 S .3 + (K + A2>s 2 + (BiK + Ai>S + BqK + Aq 

S'+ + AgS^ + A2S2 + AjS +. Aq 

The roots of this equation which are the roots of the numerator are seen to equal the 
poles of the closed-loop transfer functions T^. Because the locations of the 
closed-loop poles in the complex plane determine the stability of the dynamical 
system being considered, these will be plotted as a function of the parameter K. As 
in the open-loop case, the eigenvalues can be determined by using the root finding 
program. However, before the program can be implemented, the coefficients must be 
replaced by their numerical values. These have been introduced by the use of the 
previously described substitution routine. Subsequent to the introduction of the 
numerical values, the root finding program was used to compute a set of close-loop 
poles. The range of values obtained was sufficient to permit the plotting of the loci 
for the two dominant modes of motion. The tabulated results are shown in table 1, and 
the eigenvalue loci are plotted in figures 3 and 4. 

TABLE 1.- CLOSED-LOOP EIGENVALUES 


K 

Phugoid mode 

Short period mode 

0 

-0.00996 

± 

3 

0.1632 

-0.8403 

±j 

1.3839 

1.0 

-.0962 

± 

j 

.1457 

-.7541 

±j 

1.6710 

2.0 

-.1467 

± 

j 

.1014 

-.7035 

±j 

1.9343 

3.0 

-.1779 

± 

j 

.0258 

-.6723 

±3 

2.1722 

3.1 

-.16695 



-.1938 

-.6699 

±j 

2.1948 

3.2 

-.1510 



-.2145 

-.6675 

±3 

2.2171 

3.3 

-.1424 



-.2276 

-.6652 

±3 

2.2392 

3.4 

-.1362 



-.2382 

-.6631 

±3 

2.2611 

3.5 

-.13123 



-.2474 

-.6610 

±3 

2,2829 

3.6 

-.1271 



-.2555 

-.6589 

±3 

2.3044 

3.8 

-.1205 



-.2698 

-.6551 

±3 

2.3470 

4.0 

-.1154 



-.2821 

-.6515 

±3 

2.3888 

10.0 

-.0775 



-.4215 

-.6007 

±3 

3;4188 

20.0 

-.0686 



-.4750 

-.5785 


4.6567 

50.0 

-.0637 



-.5101 

-.5634 

±j 

7.1892 

100.0 

-.0621 



-.5224 

-.5580 


10.0839 

1,000.0 

-.0607 



-.5337 

-.5531 

±3 

31.6494 

10,000.0 

-.0605 



-.5349 

-.5526 

±j 

100.0084 

100,000.0 

-.0605 



-..5349- 

-.5525 

±j 

316.2304 


ZEROS : -0.0605; -0.535 ; ±» 
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OPEN LOOP 


u « 1.62 rad/sec 
r = 0.52 



t 4 


^3 


i2 


Hi 


1 1 1 ^0 


ROOT LOCUS FOR 
SHORT PERIOD MODE 


I I L. 

4 3 2 


Hi 



H2 


J4 

0 


Figure 3.- Root locus for short-period mode. 
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CONCLUSIONS 


It has been demonstrated that the technique of symbolic computation can be used 
to facilitate the formulation and processing of the equations of motion of dynamical 
systems. The manual formulation of these equations, which involves a variety of 
matrix operations, including Jmatrlx inversion, is time consuming and subject to human 
error. An example is included which demonstrates the facility with which the system 
dynamics matrix and the control distribution matrix from the state space formulation 
of the equations of motion can be processed to obtain eigenvalue loci as a function 
of a system parameter. The example describes a single-input, single-output state 
variable feedback system. It should be emphasized, however, that the technique is not 
limited to the study of systems of this type, but can be applied with equal facility 
to the study of more general systems. 
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